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We describe in this article a new code for evolving axisymmetric isolated systems in general 
relativity. Such systems are described by asymptotically fiat space-times which have the property 
that they admit a conformal extension. We are working directly in the extended 'conformal' manifold 
and solve numerically Friedrich's conformal field equations, which state that Einstein's equations 
hold in the physical space-time. Because of the compactness of the conformal space-time the entire 
space-time can be calculated on a finite numerical grid. We describe in detail the numerical scheme, 
especially the treatment of the axisymmetry and the boundary. 



(N 
O 
O 



<N ; 

> ■ 

OS 

o 
o 

(N . 
o ■ 

"o 

I ' 



X 



I. INTRODUCTION 

Currently, there are three main approaches to solve 
Einstein's field equations using numerical methods. They 
differ in the way they set up systems of equations and the 
appropriate initial boundary value problems which they 
try to solve. The main approach is based on variants of 
the ADM equations which are written as a hyperbolic 
system of evolution equations. The IBVP is set up on 
space-like hypersurfaces with boundary 'far out' in a re- 
gion where it is assumed that the space-time is such that 
in a certain approximation it is close to the Minkowski 
or Schwarzschild space-time or a perturbation thereof. 
Radiation is extracted by matching to such approximate 
space-times. 

The second approach is based on the characteristic ini- 
tial value problem which is based on outgoing null hyper- 
surfaces. The Einstein equations can be written in this 
context in a hierarchical form which lends itself very eas- 
ily to numerical treatment. The use of outgoing null hy- 
persurfaces has the advantage that radiation of any type 
and, in particular, gravitational radiation can be followed 
all the way out to an asymptotic observer (LIGO) located 
at null-infinity. By using a compactified radial coordinate 
it is possible to read of the radiation at finite points. 

There are some attempts to match these two ap- 
proaches in the so called Cauchy-Characteristic Match- 
ing procedure (CCM) where the two schemes exchange 
appropriate information across an interface to provide 
boundary conditions at a time-like hypersurface (see ||] 
for a recent review). 

The third approach, termed the conformal approach, 
is somewhat different. While the previous methods set 
up equations directly for the physical quantities in the 
physical space-time in this approach one formulates the 
so called 'conformal field equations' for quantities on a 
conformally related space-time, which state in an indirect 
way, that the physical Einstein tensor vanishes. This in- 
direction introduces some freedom which is used to com- 
pactify the space-time in so that it is possible to hold 
the entire (semi-global) problem on a finite grid. In this 
way one can solve a global problem with finite resources. 
Since the whole space-time is accessible one can also read 



off radition without further approximation. Due to these 
properties the conformal approach is well suited to study 
global properties of space-times. 

All these approaches have certain more or less severe 
draw-backs. Common to all of them is the fact that the 
constraints which can be satisfied initially very accurately 
are violated during the course of the evolution and this 
'constraint violation' seems to grow in all codes with an 
exponential rate. 

The standard Cauchy approach suffers from the fact 
that the outer boundary is very difficult to handle. This 
is due to a lack of understanding of the appropriate IBVP 
in the sense that there is no boundary condition known, 
which would be physically meaningful (modelling the fact 
that the space-time should be asymptotically fiat) as well 
as numerically stable. The characteristic codes crash 
when the underlying outgoing null hypersurfaces form 
caustics and self intersections. In the conformal approach 
there are no (known) difficulties out of principle. The 
problems arising there as we will point out in this work 
(see also a recent discussion in j^]) are of a numerical 
nature or due to the complexity of the equations. 

The conformal approach has been pioneered by Hiibner 
who implemented the first code treating a spherically 
symmetric space-time with a scalar field coupled to grav- 
ity He also succeeded to evolve general initial data 
which were close enough to fiat data to evolve into a 
global space-time with a regular i+ [|j (see also ||]). 

We present in this work a code to solve the conformal 
field equations for space-times which are axisymmetric 
with a regular axis. In contrast to previous work 0, |[ ||] 
where the unphysical assumption was made that there 
are no fixed points now we assume the existence of a reg- 
ular axis which implies that we can now treat physically 
relevant isolated systems. 

Our motivation for developing a 2D instead of a full 
blown 3D code is that we want to have a code with rea- 
sonable turn around times that allows us to study the 
effects of changing various properties such as gauges, 
boundary conditions etc. within a reasonable time. 
Futhermore, in contrast to a 3D code a 2D code allows 
us to compute the space-times with a much higher reso- 
lution. 

The plan of this paper is as follows. In sec. II we review 
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the basics of the conformal approach, write down the con- 
formal field equations and discuss the implementattion of 
the axisymmetry. Sec. ID is devoted to the description 
of the numerical methods used in the code. In partic- 
ular, we discuss the evolution algorithm, the boundary 
conditi ons and the issues arising from the axisymmetry. 
In sec. IV we present the results obtained so far and con- 
clude with a summary in sec. 0. 



II. THE GENERAL SETTING 
A. The Conformal Field Equations 

In this section we first give a short introduction into the 
conformal picture and provide the conformal field equa- 
tions. Then we introduce the hyperboloidal initial value 
problem for the conformal field equations by making a 
3 -|- 1-decomposition with respect to a frame into a sys- 
tem of constraint equations and a system of symmetric 
hyperbolic evolution equations. 

We then describe how to implement the axisymme- 
try with special emphasis on the regularity of the axis. 
The symmetry reduces the problem from a 3D to a 2D 
problem, but does not lead to a significant reduction of 
components. We therefore assume additionally, that the 
Killing vector field , which generates the axisymmetry 
is orthogonal to the hypersurfaces of constant 0, where 
(j> is the parameter along the orbits of (i.e., ^(0) = 1). 
This gives us an additional discrete symmetry, which re- 
duces the number of components from 65 for the general 
equations to 35 independent components. 

Suppose we have a vacuum solution (A4, g) of the Ein- 
stein equations, which we will call the 'physical space- 
time'. This space-time is asymptotically flat at null- 
infinity ('light-like asymptotic flatness') and hence de- 
scribes an isolated general relativistic system, if it admits 
the following construction. 

There exists a smooth Lorentz manifold {M,gab) and 
a diffeomorphism : — > A^, such that: 

1. ^(A4) C M and 'i'{M) has a smooth boundary, 
d^{M) = ^with J^= J^UJ^ and J^nj^ = 0. 

2. there exists a smooth function fl on Ai, such that: 
gab - nH^-^)*gab andn^O on *(Af), 

3. on the conformal factor il = and dil ^ 0, 

4. every null geodesic in A4 has a future endpoint on 

and a past endpoint on , 

5. the Ricci tensor Rab of the physical metric gab van- 
ishes in a neighborhood of Rab — 0. 

This construction states, that the physical manifold 
(M^g) is conformal to the interior of the 'unphysical' 
manifold (A4,g). In order to simplify the notation, we 



will discard further on the diffeomorphism ^' and iden- 
tify the points p G A4 and ^{p) £ M. The boundary 
^ of A4 in M. represents 'infinity' of the physical man- 
ifold. The asymptotic region ^ is in a. sense unique for 
all asymptotically flat space-times and defines a unique 
'background' structure, on which one can define the mass 
of a space-time and the gravitational radiation measured 
by an observer at infinity. 

Since one has immediate access to the geometric struc- 
ture at J^in the unphysical manifold Al, it is reasonable 
to solve the conformally transformed Einstein equations 
directly in Ai. Now the problem arises, that the trans- 
formed Einstein equations are formally singular at ^. 
This problem was solved by Friedrich [|lO| , who proposed 
the conformal field equations, a system of regular equa- 
tions for the geometric structure on . If there exists a 
triple (A^,(7,r2) satisfying the conformal field equations 
then the pair {M,g), where M — 0] and g = il^^ g 

is a vacuum space-time. 

Let us consider the abstract manifold Jv[ with the ad- 
ditional structure of a metric gab and a metric compatible 
connection Vq. T°'bc and R°'bcd are the torsion and cur- 
vature tensors, respectively, associated with Vq. Addi- 
tional variables are the conformal factor SI, a 1-form Ea, 
a scalar function 5*, a symmetric traceless tensor field 
$Qf, and a completely traceless tensor field Kabcd with 
the symmetries of the Weyl tensor. 

Then the conformal field equations read: 
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flKa 



bed 



Ta 
be 

25c[a*fc]<i - 25<i[a*6]c 



= 0, 



(1) 





-4:9c[a9b]d-^ 


-0, 


(2) 




b + 2g6[cVa]A - Kcab''''^d 


= 0, 


(3) 




'^aK'^bcd 


= 0, 


(4) 




Va^ - Sa 


= 0, 


(5) 




VaSfc - gabS + fi$a6 


^ 0, 


(6) 




h - riVaA - 2AI], 


= 0, 


(7) 




2f7S' - 2f72A - SaS" 


= 0. 


(8) 



The first equation states, that the torsion tensor T'^bc 
vanishes, so that Vq is the Levi-Civita connection of the 
metric gab and ^ is the decomposition of the curvature 
tensor into its irreducible parts. Since one can show that 
the Weyl tensor C^bcd necessarily vanishes on ^it is writ- 
ten as Cabcd = ^Kabcd- ^ab IS the traccless part of the 
Ricci tensor and A is proportional to the Ricci scalar, 
R = 24A. Eq. (I) serves to get first order equations 
and, therefore, must be regarded as the definition of E^. 
Eqs. (H) and (H) are the expressions for the conformally 
transformed traceless part of the physical Ricci tensor 
and the physical scalar curvature, respectively. Eqs. (|^) 
and (|^) are the conformally transformed Bianchi identi- 
ties of the physical manifold, which form a central part 
of the conformal field equations. 

Now we have to set up an initial value problem for 
the conformal field equations. We do this by the usual 
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3 + 1-decomposition. Since we are working in the tetrad 
formalism, we pick a tetrad with g^b = VtJ.iy^a ® 
(e^ being the dual frame of covectors) where t° = eg is 
orthogonal to the hypersm'face S of constant coordinate 
time t, so the manifold is sliced in the way 7W = R x 
S. In a tedious but straightforward procedure we get 
a system of symmetric hyperbolic evolution equations, 
together with a system of constraints. 

In this process, we also have to fix the gauge freedom. 
This is done here by adding appropriate equations to the 
system, which contain four free functions, which deter- 
mine the coordinates, six free functions corresponding to 
the freedom in the choice of the tetrad. Since the confor- 
mal field equations are invariant under conformal rescal- 

where lj has to be strictly 
positive, there is one additional function to fix the con- 
formal factor. The details of this hyperbolic reduction 
together with the fixing of the gauges are described in 
and 

Since our version of the conformal field equations is 
not suited to treat space-like infinity i*^ in an appropri- 
ate way (see for a discussion of other formulations 
of the conformal field equations) we choose as the initial 
data hypersurface a hyperboloidal slice of M.. This hy- 
perboloidal slice intersects J'^ in a space-like 2-surface, 
in contrast to an asymptotically Euclidean hypersurface 
which reaches spatial infinity i*^. Since a hyperboloidal 
slice is not a Cauchy surface for the entire space-time, 
we are not able to calculate the entire space-time, but 
we can do semi-global evolutions into the future. This 
is enough for our purposes, because we are mainly in- 
terested in the asymptotic regions near and including 
where we can read off the gravitational radiation emitted 
within the manifold. 



B. Axisymmetry 

We assume axisymmetry of our physical space-time A^, 
see for a definition. The assumption of axisymmetry 
is motivated by the fact, that we can reduce the numerical 
problem from a 3D to a 2D problem and that this sym- 
metry still allows gravitational radiation, see for the 
symmetries which admit radiation for light-like asymp- 
totically flat space-times. This allows us to simulate in- 
teresting radiating systems with high accuracy already 
on small computers. The symmetry will be transferred 
to the conformal space-time M. , if the conformal factor 
is also invariant under the symmetry C(^VL = 0, which we 
assume explicitly by partly using up the gauge freedom 
in r2. 

The reduction from a 3D to a 2D problem is achieved 
by evolving only within one of the hypersurfaces of con- 
stant (j). All these hypersurfaces are isometric and so it is 
enough to calculate all fields on one of them. The tricky 
part of axisymmetric systems is the analytic and numer- 
ical treatment of the axis. From the analytical side, one 
would like to exclude pathological things, like a massive 



string on the axis. This means we would like to have a 
space-time, which is regular on the axis. The regularity 
of the axis is usually ensured by the condition of elemen- 
tary flatness jlj] on the axis. But this condition is very 
hard to implement numerically. Fortunately, there is an 
equivalent condition of Wilson and Clarke in , which 
simply states, that if the Killing vector has the same 
form as in Cartesian coordinates and the metric is differ- 
entiable in a neighborhood of the axis, which we have to 
assume anyway for the numerics, then the axis is regular. 
This is the reason, why we use Cartesian coordinates in 
our code. We evolve in the hypersurface, where ^ = 0, 
because there the Cartesian (x, z)-coordinates coincide 
with the (r, z)-coordinates of polar coordinates. In the 
following, we will call this hypersurface the 'evolution 
hypersurface'. 

Apart from the axisymmetry we have an additional 
discrete symmetry because the Killing vector is or- 
thogonal to hypersurfaces of constant 0, so the reflection 
i across the evolution hypersurface given by i— > —cj) 
and hence i— > — is an isometry. Thus, by adapting 
the tetrad e° to the Killing vector, i.e., Cg ~ in the 
evolution hypersurface we can achieve a reduction of the 
independent tensor components. For this reduction we 
have to distinguish between tensors fields, which are in- 
variant under i or which are mapped onto their negative. 
In the first (second) class all tensor components obtained 
by contracting with an odd (even) number of will be 
zero. This reduces the problem from 65 free components 
to 35 components. 



III. NUMERICS 

In this section we want to describe the numerical 
scheme, which we used to solve the conformal field equa- 
tions. Besides the conformal field equations we tested 
the same scheme also on the scalar wave equation in the 
Minkowski space-time. The numerical grid is uniform 
and rectangular, where the outermost points are used for 
the boundary. Since we use a non-staggered grid, x = 
is one of the boundaries of our grid. All of these points 
lie on the axis of the space-time, which is not a bound- 
ary of the space-time. Therefore, we have to treat these 
points in a different way, which will be described in the 
subsection on axisymmetry and the cartoon-method. 



A. Evolution 

For the evolution the Method of Lines is used. This 
is a semi-discrete method for solving partial differential 
equations. The PDEs are discretized in space to obtain 
a system of ODEs which are then integrated in time 
along the lines in space-time given by fixing the grid 
points. We use two standard methods for solving the en- 
suing ODEs: a Runge-Kutta method and the multilevel 
Adams-Bashforth scheme, both of fourth order in time. 
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In addition, we have used the iterative Crank-Nicholson 



(ICN) scheme |17 , which is an exphcit form of the im- 
phcit Crank-Nicholson scheme. These three methods are 
all explicit, which makes them easy to implement in the 
code. Nevertheless implicit methods would have better 
stability properties, but compared to explicit methods 
they would be much slower, because one has to solve a 
big system of nonlinear, coupled equations at each time 
step. 

The spatial discretization can be chosen between sec- 
ond order and fourth order central stencils, but we use 
fourth order discretization in all test runs, as there was 
no difference in the stability properties of the method, 
only in the accuracy of the solution obtained. 



B. Axisymmetry and the Cartoon-Method 

We stated in the last chapter, that we use Cartesian 
coordinates, so that the axis is regular. Despite the fact 
that we are using non adapted coordinates, we can reduce 
the problem from 3D to 2D by evolving in the hypersur- 
face, where = and the Cartesian coordinates {x, z) 
coincide with the cylindrical coordinates (r, z) . The prob- 
lem, which arises now, is, that we will need derivatives 
across this evolution hypersurface. We solve this by ap- 
plying an idea of Alcubierre et al. |Q, which they call 
the cartoon- method. They have proposed, that instead of 
using the symmetry in adapted coordinates d(f,T — 0, one 
can use the symmetry as a tensor equation $*T = T in 
Cartesian coordinates, where T is an arbitrary axisym- 
metric tensor. Because of the use of Cartesian coordi- 
nates, we get in this way a regular axis and at the same 
time we avoid the coordinate singularity of adapted coor- 
dinates and the resulting singular components at the axis. 
One could certainly cure the singularities of adapted co- 
ordinates by a regularizing procedure for the components 
and the equations at the axis. But it should be clear, 
that this change of the equations is a possible source of 
instabilities and for our large system of 35 independent 
components, it would have been at least tedious to check 
the behavior of each component at the axis (sec , ||2^ 
for such regularization procedures and plf for a more re- 
cent approach). As Cartesian coordinates are often used 
in 3D codes, the cartoon method has the additional ad- 
vantage that we can easily compare the results of our 
code with that of 3D codes. In the following we give the 
formulation of the cartoon-method in the tetrad formal- 
ism. 

Since we are using a tetrad basis we not only have to 
ensure the regularity of our coordinates, but also the reg- 
ularity of our tetrad at the axis. Additionally, we have 
the constraint, that one basis vector e° has to be orthog- 
onal to the evolution hypersurface, so that we can use 
the discrete symmetry i to reduce the number of inde- 
pendent components. Apart from these two conditions 
we can arbitrarily fix the tetrad in our space-time. We 
first introduce an invariant tetrad where ~ and 



$*/^ = /^. Similar to adapted coordinates, this adapted 
tetrad is singular on the axis. But we have ensured, that 
our manifold is elementary flat at the axis, so we only 
have to rotate the tetrad in the neighborhood of the axis, 
to a tetrad which has a smooth limit while approaching 
the axis. So we make the ansatz, that the tetrad is given 
by the rotation = R^'^ f^, where 



Rn 



/ 1 \ 
10 


- 

p / 



^ 



Now we calculate the frame = /^^S" with respect to 
the frame on the evolution hypersurface, where 9" are 
the coordinate vectors: 

We get the following relation between the invariant tetrad 
/° at a point p on the evolution hypersurface and a point 
^{p) on the orbit of the isometry through p: 



(9) 



And together with the ansatz of our tetrad ~ Rf^'^ f^, 
we get the tetrad components anywhere in the manifold 
with respect to the components on the evolution hyper- 
surface: 



c,^<i>{p))^R.^mp))c/{p)<i>, 



(10) 



Now for regularity of the tetrad, we only have to ensure 
that there exists a uniform limit for the c^'^, when p 
approaches the axis. The coordinate vectors themselves 
have a regular limit because of the usage of Cartesian 
coordinates. The necessary behavior of the components 
c^'^ at the axis is then built implicitly into the code. 

The same procedure applied to an arbitrary tensor 
T°''' "cd... with respect to the tetrad basis results in the 
following relation for the components: 

{<p*Tr'^^-mp)) = iR-'),rmp)) ■ . .r^^^^-(p).(ii) 

All occurring geometric objects T are axisymmetric 
= T, so we get the following relation: 



{x, y, z) = (i?-i)^, {x, y,z)... r^^^-( V^^T^, 0, zI12) 



In order to ensure a smooth limit for the tensors, this 
formula leads to conditions for the components on the 
axis, which are also built implicitly into the code. 

This formula is also used to calculate the derivative 
which points away from the grid. The problem here 
is, that for the derivative at grid point {xi,Zj) of the 
evolution hypersurface in y-direction, we need the ten- 
sor components at points (\/ xf + ■m'^{Sy)^, 0, zj), where 
m depends on the employed discretization. These points 
will in general not coincide with a grid point. So one has 
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to use interpolation to calculate the components there. 
In the paper on the cartoon-method by Alcubierre et al. 
psf they propose to use Lagrange interpolation. How- 
ever, in 1^ wc show for the scalar wave equation, that 
this interpolation method results in an unstable evolution 
scheme and propose alternative schemes, which are sta- 
ble. Unfortunately even with these alternative schemes 
we encounter instabilities in the case of the conformal 
field equations. 



C. Boundary 

In our current approach, we set the boundary of our 
numerical grid into the unphysical part of the conformal 
manifold. As the physical boundary c^is a characteristic 
hypersurface, analytically no influence from the unphys- 
ical part can enter the physical part. Our results show, 
that this is also true in our numerical treatment. So 
we do not have to provide physically correct boundary 
conditions, like e.g. constraint preserving boundary con- 
ditions, as one should if the boundary would be in the 
physical part. Instead we only have to impose numerical 
boundary conditions, which lead to a stable and conver- 
gent numerical scheme. 

This is not a much simpler task. We follow a standard 



characteristic fields of the following form: 



^ 3 ^^out 



scheme |22 to get numerical boundary conditions for a 



symmetric hyperbolic system, which imposes boundary 
conditions on the ingoing local characteristic fields. The 
outgoing local characteristic fields are determined by the 
data on the grid and are extrapolated to the boundary 
point. We only provide second order boundary condi- 
tions, which are simpler and less computationally expen- 
sive, than fourth order conditions. This has the obvious 
disadvantage, that we can use only a second order spatial 
discretization at the last inner grid point, which results 
in a global third order scheme, as the numerical results 
indicate. 

In order to get the local characteristic fields, one con- 
siders only the principal part of the evolution equations 
t°'daU — A^efdaU. Then one determines the outward 
pointing unit normal vector n° to the boundary and de- 
composes the principal part into the parallel and or- 
thogonal part to the normal vector: dtu = {A^rii — 
t°-naE)n°-daU + tangential part, where E is the identity 
matrix. The tangential part can be discarded, because it 
is only important, what moves across the boundary. We 
determine the eigenvalues A and eigenvectors 11 of the 
matrix A^rii — f^UaE and classify the local characteris- 
tic fields n, by the sign of the corresponding eigenvalue 
(this is the reason, why we have to ensure, that the nor- 
mal vector is outward pointing). 



A > 
A < 
A = 



mgomg, 
outgoing, 

tangential to the boundary. 



Now we impose boundary conditions on the ingoing local 



where the matrix K^j has to fulfill < 1 and the 

function g^{t) can be specified with certain restrictions. 
The values of the outgoing characteristic fields are ex- 
trapolated from the data on the grid. The derivation 
of the boundary conditions seems to indicate, that we 
have to extrapolate in the direction of the normal vector. 
This is in general not well defined in a manifold and as 
the normal vector has in general also time components, 
it may even require to extrapolate in time not only in 
space, which is obviously impossible. Therefore we take 
the most simple approach and extrapolate along the grid 
lines. 

The matrix R^j acts like a reflection coefficient. As 
the boundary has no physical meaning and is simply the 
boundary of the numerical grid, reflection of the outgoing 
characteristic fields would be wrong and we use highly 
absorbing boundary conditions, by setting R^j = 0. 

At the moment we are not able to generate general 
initial data. For this reason and because we want to 
test our numerical scheme we use initial data generated 
from exact solutions. Therefore, we set to the value 
computed from the exact solution, when it is not stated 
otherwise. 

So far this has been the standard procedure for pro- 
viding boundary conditions, if a fixed background metric 
is available. In our case, the metric is determined by the 
tetrad components, for which we also have to provide 
boundary values. This means in principle we cannot de- 
termine the normal vector at the boundary, because we 
do not know the metric at the boundary. We solve this 
problem by extrapolating the necessary tetrad compo- 
nents. From the extrapolated tetrad components we cal- 
culate the eigenvalues of the characteristic fields. This 
should not result in a significant error, because only the 
sign of the eigenvalue is important. A bit more tricky is 
the transformation from the characteristic fields to the 
physical fields. This yields nonlinear coupled equations 
for the tetrad components, which occur in the normal 
vector. We solve the occurring equations with the New- 
ton method. 

The scheme, we have presented for providing bound- 
ary conditions, is a generalization of the standard scheme 
of | |23[ |, where one allows varying coefficients in the sym- 
metric hyperbolic principal part and in addition works 
without a background metric. Therefore it is not clear, 
whether the theorems derived for the standard case are 
still valid here. 



IV. RESULTS 

In this section we describe the results, which we have 
obtained so far. At first we describe the results for the 
scalar wave equation, before we examine the conformal 
field equations. In both cases we use the same numerical 
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scheme. This means only the equations are changed, the 
evolution and the treatment of the boundary and the axis 
are the same. 

In all cases the grid size s is: s — Nx{2Nx — 1), where 
Nx is the number of grid points in x-direction, and Nz — 
2Nx — 1 the corresponding number of points in the z- 
direction. This ensures, that the origin lies on a grid 
point. 

Errors Q are measured with the Li-Norm: 

X! X! \\Hx{n),z{m),t) - Un^„i{t)\\i, 

n—l m—1 

where u is the exact value of u at {x{n), z(m)) and Un,m 
the value at the grid point {n,m). In these sums the 
boundary points are excluded. Here u stands for either 
the vector of all variables or the vector of all constraint 
components. 

A. Axisymmetric scalar wave equation 

At first we would like to discuss shortly the results 
of our code for the axisymmetric scalar wave equation 
in flat space-time. As we have written the scalar wave 
equation also in symmetric hyperbolic form, we can use 
it as a testbed for our numerical scheme, especially for 
the treatment of the boundary and the axisymmetry. 

For the scalar wave equation Lagrange interpolation 
leads to an unstable numerical scheme. In |3| we in- 
troduce a second order interpolation procedure for the 
cartoon method, which results in a stable second order 
evolution scheme. We have also found a fourth order 
interpolation scheme, which seems to be stable, but we 
cannot prove it yet. The resulting evolution scheme is 
third order accurate. This is a result of our second or- 
der boundary treatment, which changes our fourth order 
scheme globally to a third order scheme. 

The code was also tested on robust stability, which was 
proposed by Winicour. This tests the scheme by evolving 
random initial data with random boundary data whereby 
all possible solutions of the numerical scheme will be ex- 
cited. In our numerical tests we see no emerging devia- 
tion from the mean after one million timesteps, where we 
stopped the evolution. 

These results show, that our numerical scheme is a suit- 
able method for solving axisymmetric systems as long as 
one can find an interpolation method, which is compati- 
ble with the evolution scheme. 



B. Axisymmetric Conformal Field Equations 

Unfortunately the very good results from the scalar 
wave equation do not transfer to the conformal field equa- 
tions. It was not possible to get a stable evolution for any 
of the space-times, we have tested. But this does not 
mean, that our results are physically meaningless. For 



each space-time, where a regular point i"*" exists, it was 
possible to evolve until i"*" and sometimes much further 
on. 

As we have no general initial data yet, we tested our 
code with exact solutions. The exact expressions for the 
fields were calculated by GRTensorll and then trans- 
lated directly into C code. 

1. Minkowski Space-Time 

At first it may seem trivial to test the code with the 
Minkowski space-time. But the conformal compactifica- 
tion of Minkowski space-time has curvature and can even 
be time dependent, if one uses a time dependent confor- 
mal gauge. We have tested three different compactifica- 
tions of the Minkowski space-time. 

The first one is the standard compactification, as it is 
introduced in textbooks: 

— cos(i) -I- cos(r), 
gab = dt^ — dr^ ~ siv? {r)duj^ , 

where dw^ is the line element of the unit sphere. The sec- 
ond one is a rescaling of the standard compactification, 
which has the property, that is not a finite point: 

n = cos(i)(l + r2) + l-r2, 

gab — —{l + r^)dt'^ — dx^ — dy^ — dz'^, 

and the last one is the static compactification 

gab — ^{r^ - l)'^dt'^ + (xdx + ydy + zdz)dt 
~dx^ — dy^ ~ dz^ . 

a. Convergence We show convergence of our numer- 
ical scheme for the standard compactification. In Fig. |l| 
we plot the error summed over all fields and the con- 
straint violation for the whole grid and for the physical 
part, O > 0, at coordinate time t — 2.1 for increasing 
number of grid points N, the initial time is tmi = 1.8. 
This yields a convergence order for the total error of 
2.903 ± 0.002, which is in good agreement with our ex- 
pectation of a third order scheme. The total constraint 
violation converges to zero with order 3.32 ± 0.08. The 
physical errors are both much smaller than the total er- 
ror, which comes from the fact, that the error is concen- 
trated at the boundary, which does not lie in the physical 
space-time. 

A problem we face at the moment is, that the be- 
haviour of the error does not generalize to higher grid 
dimensions. There fast growing instabilities arise at the 
boundary, which are not present at lower resolutions. 
These instabilities dominate the total error after a short 
time and lead to a crash at the boundary. 
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is stable under small perturbations entering from the 
unphysical part. Second the noise from the boundary, 
which is also added to the Weyl tensor, enters the grid, 
but does not enter the physical space-time. Note that 
the scale of the plot is 10^^ and the plot is cut off at 
1.2 X 10"^. The plot would have been 1000 times bigger 
to cover the whole range of C. So is indeed a charac- 
teristic hypersurface not only analytically but also with 
great accuracy numerically. Near at t — tt a, small 
amount of noise enters the physical space-time, but this 
seems to be related to the very poor resolution of the 
physical space-time in this area. 

Further support for the property of ^ can be gained 
by the following Fig. || of the error and constraint vio- 
lation for the whole grid and the physical part. It can 



FIG. 1: Total error with respect to the exact solution for the 
standard compactification, for A^=20,40,60,80 and 100 points 
in the x-direction 

b. Numerical Properties of The great advantage 
of the conformal approach, that the boundary of the grid 
does not lie in the physical space-time, in opposite to the 
other approaches in numerical relativity, can be further 
illustrated with the following example. We start with 
initial data from the standard compactification and add 
during the evolution uniform random noise to the exact 
values for the ingoing characteristic fields at the bound- 
ary. The maximum of the noise is (5 = 0.005, we use 
40 X 79 grid points, initial time is tini = 1.8. We plot 
in Fig. I the quantity C = \EabE''^ ^ BabB"''] composed 
out of the electric and magnetic part of the Weyl ten- 
sor. In an undisturbed evolution C remains zero all the 
time, because Eat = Bab = in the initial data and this 
is preserved through the numerical evolution. In Fig. 




FIG. 2: C = \EabE''^ + BabB"-^] for the standard compacti- 
fication with noise entering through the boundary. The blaclt 
line indicates ,f (Q. = Q). 

one can see two crucial points. First is not affected 
at all by the noise from the boundary. This means, that 
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FIG. 3: Error and constraint violation for the standard com- 
pactification with noise entering through the boundary. 

clearly be seen, that the physical error is much smaller 
than the total error, the same is valid, but not with the 
same accuracy, for the constraint violation. The reason 
for this might be, that one calculates the physical con- 
straint violation at points, where f2 > 0, but through the 
discretization for the derivatives one gets an influence 
from points outside of the physical space-time. Both the 
physical error and the constraint violation in the physical 
part grow rapidly, while approaching «+ at t = tt, which 
has the reason in the low number of physical grid points 
shortly before i^ . This behavior of the error and the con- 
straint violation is not limited to this special case. Also 
for all other tested space-times, one can see, that the 
error in the physical part is small relative to the error 
in the unphysical part. We have made several movies, 
where one can see, that the error stays outside of ^. 
These movies are available on our home page ||2^ . 

c. Reproducing the static Einstein universe In Fig. ^ 
we plot the conformal factor of the standard compacti- 
fication, where we give the exact values for the ingoing 
characteristic fields during the evolution. It can be seen, 
that the physical space-time bounded by shrinks to a 
point so it is no problem to reach i+ and continue in 
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FIG. 4: The conformal factor Q for the standard compacti- 
fication at 2: = 0. The black hne indicates ^ {Q — 0) 



also very peeked solutions, which leads to a more rapid 
exponential growth of the error. 

These results are illustrated in the next plot (Fig. ||). 
As the exponential growth is different for the chosen 
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a smooth way. In principle it is now physically meaning- 
less to continue the evolution, because the whole physical 
part has been calculated. But from the numerical point 
of view it is interesting to see, for how long it is possi- 
ble to evolve further, before the code crashes. The plot 
shows, that it is possible to calculate additional one and 
a half copies of Minkowski space-time. As the Minkowski 
space-time is embedded in the static Einstein universe, 
we obtain nothing else but a part of the static Einstein 
universe. 

d. Conformal Gauges As we have pointed out in II 
in the conformal picture there is one additional gauge 
freedom apart from the coordinate gauge and the tetrad 
gauge. We are allowed to rescale the conformal factor and 
the metric according to fl ^ ujfl,gab ^'^dab, where uj 
has to be strictly positive, so that the asymptotic struc- 
ture is conserved by this conformal rescaling. 

For the second conformal manifold ([l3| ) we tested dif- 
ferent conformal gauges. Some of these gauges lead to a 
very fast exponential growth of the error and were there- 
fore useless. The most successful choice was the gauge: 

1 



1 + + z'' 



which made many components spatially constant. It is 
therefore not too surprising, that this solution had the 
smallest exponential growth of the error. Motivated by 
this gauge, we tried another time dependent one, which 
should model a nearly constant factor lu in the physical 
space-time and then a decaying factor in the unphysi- 
cal space-time in order to have very small values at the 
boundary: 

1 



FIG. 5: The effect of different conformal gauges on the total 
error for the second compactification of Minkowski ( p^ space- 
time, grid: 40 x 79, Uni = 1.3 

gauges, it seems, as if different exponentially growing so- 
lutions are excited. We do not know, whether this differ- 
ence emerges already on the level of the conformal field 
equations or on the level of the numerical implementa- 
tion. As the exponential growth is a serious problem for 
long term runs, this has to be further examined in the 
future. 

e. Comparison of the three compactifications Before 
we conclude the section on the Minkowski space-time, we 
would like to illustrate the long term behavior of the three 
different compactifications in Fig. The absolute value 




This function has the serious drawback, that it forms a 
sharp peek at the axis as time increases. This results in 



FIG. 6: The long term behavior of the total error of three dif- 
ferent compactifications of Minkowski space-time, the marks 
are set every 1000 time steps, grid size: 40 x 79 

is not comparable, because the step size is not the same 
for the three cases. But the exponential growth can be 
compared, because in principle this should not depend on 
the step size. As expected, the fully static solution runs 
for the longest time. The big difference in the stability 
between the standard compactification and the light-like 
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compact one seems to arise, because of the greater spatial 
variance of the hght-hke compact solution. 



2. Boost- axisymmetric solution of Bicdk and Schmidt 

As a further testbed for our code, we chose an 
asymptotically flat solution from the general class of 
boost-axisymmetric space-times analysed by Bicak and 
Schmidt |2q| . Since we are working in the conformal pic- 
ture it is necessary, that the solution admits a regular ^ 
and is regular on the axis. Indeed there is a whole fam- 
ily of such solutions in the class of boost-axisymmetric 
solutions, which fulfills these conditions, except for two 
singular points on , see This family of solutions 
describes two self-accelerating particles and the singular 
points on are the places, where the particles hit . 



{xdx + ydyY {tdt — zdzY 



+e 



ydxY 



^^(tdz_ 



zdtf 



<2 



-t^ + + + z^ 



with 



V. SUMMARY 

In this article we have presented a new code for simu- 
lating axisymmetric, isolated systems in general relativ- 
ity. We have implemented the axisymmetry by adapting 
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FIG. 7: Error and constraint violation for the boost- 
axisymmetric solution with strong data m = 1.75, grid size: 
40 X 79 



- z^) = m{t^ - z'^ - - y'^), 
X{x^ +y^,t^ - z^) = m^{x^ -\-y^ f ~m{x^ + y^ + - z^), 

where to is a free parameter. In the limit m — > the 
solution approaches Minkowski space-time. In order to 
avoid the singularities in the evolution we chose our hy- 
perboloidal initial data hypersurface in the future of the 
singular points and so we are able to calculate the rest of 
the space-time up to i'^ . 

We first show an error plot in Fig. |^ for a run with 
TO = 1.75 and a grid of 40 x 79 points, which shows, that 
there is no problem to reach i'^ a,t t — Q and to evolve 
further. The high difference between total and physical 
error arises, because of the very high values near the grid 
boundary. This error does only enter the physical space- 
time to a certain degree, which shows again the advantage 
of the conformal approach to have the boundary in the 
unphysical part of the space-time together with a charac- 
teristic hypersurface J^, which separates the unphysical 
from the physical part. 

To show that this case is non-trivial in the sense that 
the curvature does not vanish we plot in Fig. ^ the quan- 
tity C = \EabE''^ + BabB'"'']. The scale of the plot in- 
dicates, that we are really dealing with data in the non- 
linear regime and that the numerical scheme can deal 
with large gradients in the data. 
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the cartoon method to the tetrad formalism and formu- 
lated boundary conditions, which are compatible with 
the symmetric hyperbolic principal part of our evolution 
equations. We have tested the numerical scheme by ap- 
plying it to the axisymmetric scalar wave equation. The 
results show, that we can solve the wave equation in a 
stable and convergent way, by choosing a suitable inter- 
polation scheme for the cartoon method. The fact that 
the code is still unstable for the quasilinear conformal 
field equations needs to be examined in further detail. 

The numerical test for the conformal field equations 
with hypcrboloidal initial data from different compact- 
ifications of the Minkowski space-time and strong data 
from a family of boost-axisymmetric space-times show, 
that the semi-global simulation of space-times until i"*" 
and further on is possible. In addition, we have demon- 
strated, that ^ is also numerically a characteristic hy- 
persurface and acts as a boundary between physical and 
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unphysical space-time with good accuracy. 

There remain several tasks. We need to develop mod- 
ules to extract the radiation on ^ and to determine the 
Bondi mass. These are straightforward implementations 
of the algebraic transformation between our tetrad and 
the tetrad adapted to Furthermore, we need to im- 
prove our boundary conditions to be at least third order 
accurate so that we do not loose one order of global ac- 
curacy as we do now. An analysis of the initial boundary 
value problem needs to be carried out (possibly along 
the lines of pi]) in order to be able to provide bound- 



ary conditions which are compatible with the constraints. 
This would minimize the influence of constraint violating 
modes generated at the boundary. 
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